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Abstract 

The  lattice  contribution  to  thermal  conductivity  of  single-walled  carbon  nanotubes  with  three  different  screw  symmetry  (chirality)  is  studied 
using  the  Green-Kubo  relation  from  linear  response  theory  and  molecular  dynamics  based  thermal  current  auto-correlation  functions.  The 
interactions  between  carbon  atoms  are  analyzed  using  the  Adaptive  Intermolecular  Reactive  Empirical  Bond  Order  (AIREBO)  potential.  The 
results  obtained  show  that,  due  to  an  exponential-decay  character  of  the  long-time  thermal  current  auto-correlation  functions,  quite  accurate 
lattice  thermal  conductivities  can  be  obtained  using  computational  cells  considerably  smaller  than  the  phonon  mean  free  path.  In  addition,  the 
computed  lattice  contributions  to  thermal  conductivities  are  found  to  agree  within  a  factor  of  two  with  their  counterparts  obtained  using  the 
Boltzmann  transport  equation.  Also,  chirality  is  found  to  affect  lattice  thermal  conductivity  by  as  much  as  20%. 
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1.  Introduction 

As  the  size  of  electronic  and  mechanical  devices  is  re¬ 
duced  into  nanometer  length  scale,  materials  thermal  con¬ 
ductivity  becomes  an  increasingly  important  since  the  oper¬ 
ations  of  such  devices  may  typically  require  that  significant 
amounts  of  heat  be  dissipated  in  a  small  region.  In  general, 
however,  experiment  determination  of  thermal  conductivity 
in  nano-scale  devices  is  quite  difficult  (e.g.  [1]),  particular 
in  the  case  of  devices  with  complex  geometries.  Fortunately, 
reliable  theoretical  and  computational  methods  have  been 
developed  for  predicting  the  thermal  properties  of  nanoscale 
materials  and  devices  (e.g.  [2]). 

Because  of  a  remarkable  combination  of  their  properties 
(e.g.  high  hardness  and  stiffness,  light  weight,  special  elec¬ 
tronic  structure,  high  stability,  etc.),  carbon  nanotubes  are 
being  considered  as  prime  candidate  materials  for  nano-scale 
device  applications.  Consequently,  considerable  effort  has 
been  invested  in  characterizing  properties  of  carbon  nan¬ 
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otubes,  particularly  their  electronic  and  mechanical  proper¬ 
ties  [3-6].  Surprisingly,  despite  the  importance  of  thermal 
management  in  nano-scale  devices,  there  has  been  relatively 
little  progress  in  characterizing  thermal  conductivity  of  car¬ 
bon  nanotubes.  This  is  partly  due  to  challenges  associated 
with  nano-scale  experimental  measurements  mentioned 
above,  but  it  is  also  a  result  of  technological  difficulties  of 
synthesizing  high-quality,  well-ordered  nanotubes.  Conse¬ 
quently,  theoretical  calculations  of  thermal  conductivity  of 
carbon  nanotubes  are  presently  very  essential. 

Theoretical  calculations  of  thermal  conductivity  of  ma¬ 
terials  can  be  classified  as  two  main  approaches:  (a)  first 
principles  based  atomistic  simulations  (e.g.  [7-9]).  This  ap¬ 
proach  is  particularly  useful  for  nano-scale  devices  where 
the  experimental  determination  of  the  thermal  conductiv¬ 
ity  is  quite  challenging;  and  (b)  continuum  calculations 
based  on  transport  theories  such  as  the  Boltzmann  transport 
equation  [10-12].  The  main  advantage  of  the  continuum 
approach  is  that  it  enables  an  analysis  of  relatively  large 
systems.  However,  the  approach  entails  the  knowledge  of 
certain  parameters  such  as  phonon  relaxation  time  and 
phonon  density  of  states  which  must  be  determined  using 
either  experimental  measurements  (may  be  difficult  in  the 
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case  of  nano-scale  devices)  or  by  theoretical  calculations. 
An  addition  shortcoming  of  the  continuum  approach  is  that 
solving  the  governing  integro-differential  transport  equation 
may  be  quite  difficult  in  some  cases. 

Because  of  the  aforementioned  limitations  of  the  contin¬ 
uum  approach,  the  first  principles  based  atomistic  simula¬ 
tions  are  increasingly  getting  more  attention  as  a  means  of 
predicting  thermal  properties.  Besides  not  requiring  the  prior 
knowledge  of  the  model  parameters,  atomic-scale  calcula¬ 
tions  enable  quantification  of  the  effect  of  microstructure 
(e.g.  phase  interfaces  and  surface  reconstruction)  on  ther¬ 
mal  properties.  Furthermore,  atomistic  simulations  can  be 
used  to  determine  the  parameters  for  the  continuum  mod¬ 
els  discussed  above  and,  thus,  help  bridge  gap  between 
atomistic-scale  and  continuum-level  calculations. 

There  are  generally  two  main  components  of  thermal  con¬ 
ductivity  of  a  material:  (a)  an  electronic  component  which 
is  controlled  by  the  electronic  band  structure,  electron  scat¬ 
tering,  and  electron-phonon  (lattice  vibrations)  interaction; 
and  (b)  a  lattice  component  which  is  mainly  controlled  by 
phonons  and  phonon  scattering.  In  the  present  paper,  only  the 
lattice  contribution  to  thermal  conductivity  of  carbon  nan¬ 
otubes  is  considered  using  atomistic  simulations.  It  should 
be  noted  that  the  electronic  contribution  to  thermal  conduc¬ 
tivity  is  very  small  and  can  be  neglected  only  in  materials 
with  large  band  gaps.  As  far  as  carbon  nanotubes  are  con¬ 
cerned,  the  size  of  their  band  gap  is  found  to  be  dependent 
on  their  screw  symmetry  (chirality),  as  well  as  on  their  di¬ 
ameter  and  length.  The  largest  band  gap  (on  the  order  of 
1 .5  eV)  is  found  in  small  diameter,  short  (n,  m)  nanotubes 
with  the  rollup  vectors  n  and  m  satisfying  the  condition: 
\n  —  m\  ^  3 p,  where  p  is  an  integer.  For  other  types  of  nan¬ 
otubes,  the  band  gap  is  considerably  smaller  approaching 
zero  in  the  case  of  “arm  chair”  nanotubes  in  which  n  —  m. 
For  comparison,  the  band  gap  in  a  typical  insulator  is  on 
the  order  of  5-10  eV  and  of  a  semiconductor  in  a  range  be¬ 
tween  0.5  and  2.5  eV.  Based  on  these  observations,  one  may 
conclude  that  electronic  contribution  to  the  thermal  conduc¬ 
tivity  can  be  significant  in  carbon  nanotubes,  particularly  in 
those  with  metallic  behavior,  i.e.  with  a  small  band  gap. 

When  thermal  conductivity  of  solid  crystalline  materials 
is  being  calculated  using  atomistic  simulations,  due  to  a  rel¬ 
atively  large  phonon  mean-free  path  in  such  materials,  one 
might  expect  that  the  size  of  the  simulation  cell  on  the  or¬ 
der  of  several  hundred  nanometers  must  be  used  to  prevent 
phonon  scattering  from  the  cell  boundaries.  In  addition, 
quantification  of  phonon-phonon  interactions  (responsible 
for  the  finite  values  of  thermal  conductivity)  is  generally 
very  complicated  (e.g.  [13]).  In  a  recent  study,  Che  et  al. 
[2]  carried  out  molecular  dynamics  simulations  and  used 
the  linear-response  theory  based  Green-Kubo  equation  and 
energy-current  auto-correlation  functions  to  determine  ther¬ 
mal  conductivity  of  diamond.  They  found  that  while  the 
accuracy  of  the  thermal  conductivity  is  indeed  dependent 
on  the  size  of  the  periodic  cell,  an  accurate  thermal  con¬ 
ductivity  can  be  obtained  using  periodic  cells  about  40-60 


times  smaller  than  the  actual  phonon  mean  free  path.  They 
attribute  this  finding  to  the  fact  that  the  energy-current 
auto-correlation  time  is  much  shorter  than  the  energy  re¬ 
laxation  time.  The  approach  of  Che  et  al.  [2]  is  used  in  the 
present  work  to  determine  the  lattice  contribution  to  thermal 
conductivity  of  single-walled  carbon  nanotubes  of  different 
chiralities. 

The  organization  of  the  paper  is  as  follows:  a  brief 
overview  of  the  theoretical  background  related  to  the  com¬ 
putation  of  thermal  conductivity  is  given  in  Section  2.1. 
Section  2.2  and  Appendix  A  contain  details  concerning  the 
atomistic  simulation  procedure  and  the  interatomic  poten¬ 
tial  used.  The  main  results  obtained  in  the  present  work  are 
presented  and  discussed  in  Section  3.  The  key  conclusions 
resulting  from  the  present  work  are  given  in  Section  4. 

2.  Procedure 

2.1.  Theoretical  background 

At  the  macroscopic  level,  thermal  conductivity  is  defined 
from  the  Fourier’s  law  for  steady-state  heat  conduction  as: 

Tq  =  —A  ■  SIT  (1) 

where  Jq  is  the  steady-state  heat  flux  (current),  A  the  thermal 
conductivity  second-order  tensor  and  SJT  the  temperature 
gradient. 

In  general,  the  total  energy  current,  Je ,  includes  the  con¬ 
duction  heat  current,  Jq,  and  the  diffusion  energy  current, 
p.J,  where  J  is  the  particle  current,  and  p  is  the  chemical 
potential.  Thus,  the  following  relation  between  the  energy 
current  and  the  heat  current  can  be  defined  [14]: 

-  P-J.  (2) 

In  solids  under  non-extremely  high  temperatures,  the  diffu¬ 
sion  contribution  to  the  energy  current  can  be  neglected. 

In  a  system  consisting  of  discrete  particles,  the  energy 
density,  h(r),  can  be  expressed  (in  the  classical  limit)  as 
the  site  energy  of  each  particle  and  consequently,  the  heat 
current  can  be  defined  as: 

(3) 

i 

where  the  raised  arrow  denotes  a  vector  quantity,  r  the  po¬ 
sition  vector,  t  is  time  and  subscript  i  denotes  particle  i. 

One  approach  to  the  determination  of  thermal  conductiv¬ 
ity  by  atomistic  simulations  is  to  place  the  computational  cell 
in  contact  with  two  different  reservoirs  with  temperatures 
T i  and  T2  and  to  calculate  the  heat  current  when  the  system 
reaches  the  steady  state.  However,  due  to  small  dimensions 
of  the  computational  cell  (typically  10-50  nm  edge  side), 
even  a  small  temperature  difference  of  10  K  across  the  sys¬ 
tem  gives  rise  to  a  thermal  gradient  on  the  order  of  108  K/m. 
It  is  unlikely  that  the  linear  response  theory  (i.e.  the  linear¬ 
ity  between  the  heat  flux  and  the  temperature  gradient  as 


206 


M.  Grujicic  et  al.  / Materials  Science  and  Engineering  B107  (2004)  204-216 


defined  in  Eq.  (1))  would  hold  under  such  extreme  thermal 
loading.  Moreover,  this  temperature  gradient  maybe  smaller 
than  the  thermal  fluctuations  in  the  system,  making  it  diffi¬ 
cult  to  obtain  convergence  of  the  simulation  results  within 
reasonable  simulation  times. 

Due  to  the  shortcomings  of  the  approach  discussed 
above,  the  fluctuation-dissipation  theorem  from  the  linear 
response  theory  which  provides  a  connection  between  the 
energy  dissipation  in  irreversible  processes  (heat  conduction 
in  the  present  case)  and  thermal  fluctuations  (of  the  heat 
current  in  the  present  case)  of  a  system  in  equilibrium  [15] 
is  used  in  this  work.  Within  this  approach,  the  thermal  con¬ 
ductivity  tensor  can  be  expressed  in  terms  of  heat-current 
auto-correlation  functions,  CCj(t),  [15,16]  as: 

1  C°° 

<4) 


where  is  the  Boltzmann’s  constant,  T  temperature,  V  vol¬ 
ume  and 


c5(t)  =  (Jq(t)Jq(  0)>.  (5) 

Cj  (f)  is  obtained  by  phase  (particles  positions  and  momenta) 
space  (r)  averaging  as: 


(7q(t)/q(0)) 


/ d.T exp  (— j6//)/q(r)7q(0) 
f  dr  exp  (—fiH) 


(6) 


where 

H  =  YJhi.  (7) 

i 


within  the  framework  of  molecular  dynamics  simulations, 
Eq.  (7)  is  evaluated  as: 

Cj(tfj( ,o)>  =  -jJ—J2  J(t° + (») 

^corrJVcoir 

where 

t  =  At,  0  <  Nt  <  Nmd  (9) 

t0  =  No  At,  0<N0<Nmd-  Nt  (10) 


and 


Neon  —  int 


(ID 


where  t  and  to  are,  respectively,  the  current  and  the  initial 
correlation  times.  At  is  the  simulation  time  step,  Amd  is  the 
total  number  simulation  steps.  No  and  Nt  are  integers,  and 
int  denotes  operator  for  conversion  of  a  real  number  to  the 
closest  smaller  integer. 

It  should  be  noted  that  the  analysis  presented  above  is 
classical,  i.e.  no  quantum  effects  are  considered.  In  gen¬ 
eral,  quantum  effects  are  not  critical  when  the  temperature 
is  significantly  higher  than  the  Debye  temperature.  While 
no  information  is  available  regarding  the  Debye  tempera¬ 
ture  of  carbon  nanotubes,  this  temperature  is  expected  to 


be  considerably  higher  than  the  room  temperature  consider¬ 
ing  magnitudes  of  the  Debye  temperature  in  other  forms  of 
carbon.  Using  a  quantum-physics  based  analysis,  Che  et  al. 
[2]  showed  that  in  (hypothetical)  purely-harmonic  systems 
in  which  different  phonon  modes  do  not  interact,  quantum 
effects  are  negligible.  In  real  systems,  on  the  other  hand, 
phonons  are  coupled  and  this  anharmonicity  is,  in  fact,  re¬ 
sponsible  for  a  finite  phonon  mean  free  path,  and  hence,  for 
finite  thermal  conductivity.  Fortunately,  thermal  conductiv¬ 
ity  in  such  systems  is  dominated  by  low-frequency  phonon 
modes,  which  are  nearly  classical.  Hence,  a  lack  of  inclusion 
of  the  quantum  corrections  is  generally  considered  not  criti¬ 
cal  when  calculating  thermal  conductivity.  As  far  the  anhar¬ 
monicity  effects  are  concerned,  they  are  generally  accounted 
for  when  (classical)  molecular  dynamics  simulations  based 
on  realistic  interatomic  potentials  are  carried  out  and,  hence, 
such  simulations  can  be  used,  as  is  done  in  the  present  work, 
to  determine  thermal  conductivity  quite  accurately. 

2.2.  Simulation  procedure 

Molecular  dynamics  simulations  are  conducted  using 
computational  cells  which  are  of  a  finite  size  in  one  and 
infinitely  long  in  the  other  two  directions.  The  periodic 
boundary  conditions  are  applied  in  the  finite  direction. 
Each  cell  contains  a  single  ( n ,  m)  nanotube  with  the  nan¬ 
otube  axis  aligned  with  the  direction  in  which  the  cell  is 
finite.  The  dimension  of  the  cell  in  the  nanotube  direction 
is  varied  in  order  to  accommodate  nanotubes  of  different 
chirality  and,  also,  to  explore  the  effect  of  the  cell  size 
on  thermal  conductivity.  Three  types  of  carbon  nanotubes 
are  studied  in  present  work:  (a)  a  (10,  10)  armchair  type 
nanotube;  (b)  a  (18,  0)  zig-zag  type  nanotube;  and  (c)  a 
(14,  6)  nanotube.  The  first  nanotube  is  selected  because  it 
is  frequently  found  in  synthesized  nanotube  bundles,  while 
the  other  two  are  selected  on  the  basis  that  they  have  the 
diameter  compatible  to  that  of  the  (10,  10)  nanotube.  Other 
structural  characteristics  of  the  three  nanotubes  are  given  in 
Table  1  while  the  corresponding  atomic  arrangements  are 
shown  in  Fig.  l(a)-(c). 

The  interactions  between  carbon  atoms  have  been  mod¬ 
eled  using  the  Adaptive  Intermolecular  Reactive  Empirical 
Bond-Order  (AIREBO)  potential  developed  by  Stuart  et  al. 
[17].  This  potential  is  an  extension  of  the  original  Brenner’s 
Reactive  Empirical  Bond-Order  (REBO)  potential  [18]  and 
includes  non-bonding  (intermolecular)  atomic  interactions. 


Table  1 

Structural  characteristics  of  the  carbon  nanotubes  studied  in  the  preset 
work 


Nanotube 
type  (n,  m ) 

Nanotube 
radius  (nm) 

Unit  cell 
length  (nm) 

Number  of  atoms 
per  unit  cell 

(10,  10) 

1.351 

0.2477 

40 

(18,  0) 

1.404 

0.4290 

72 

(14,  6) 

1.387 

3.8130 

632 
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Fig.  1.  Atomic  configuration  associated  with:  (a)  a  (10,  10)  arm-chair  nanotube;  (b)  a  (18,  0)  zig-zag  nanotube;  and  (c)  a  (14,  6)  nanotube  of  a  general 
chirality. 


It  should  be  noted  that  while  frequently  the  interactions 
between  carbon  atoms  within  a  single  single-walled  car¬ 
bon  nanotubes  are  modeled  using  only  the  bonding  part 
of  an  inter-atomic  potential,  the  AIREBO  potential  used  in 
the  present  work  has  been  optimized  to  also  include  the 
non-bonding  interactions  between  carbon  atoms  of  the  same 
single-walled  carbon  nanotube.  Based  on  the  AIREBO  po¬ 
tential  (reviewed  briefly  in  Appendix  A),  the  total  potential 
energy  be  formally  written  as  a  summation  of  pair-wise  in¬ 
teractions  as  [18]: 

y«  =  EE^-  (12> 

i  j>i 


Consequently,  the  energy  of  site  i,  hi,  can  be  defined  as: 


However,  since  Vjj  is,  in  fact,  a  many-body  potential,  the 
equation  for  heat  current  Jq(t)  based  on  Eq.  (3)  is  somewhat 
complicated  and,  under  the  condition  of  a  zero  net  momen¬ 
tum  for  the  system,  is  given  by: 

\{t)  =  +  iE  E  •  «.•>  <14) 

i  i,j  k,l 


where 


~pkl  dVkl 

ij  dnj ' 


05) 


where  Vj  and  ,f;y  are,  respectively,  the  /-site  velocity  and  the 
i-j  sites  relative  position  vector. 

All  molecular  dynamics  simulations  are  carried  out  in  the 
present  work  using  a  fixed  1  fs  time  increment.  For  each 
simulation  run,  the  system  is  equilibrated  under  isothermal 
(T  =  300  K)  conditions  using  Berendsen  thermostat  [19]  for 
40 ps.  Subsequently,  constant  energy  molecular  dynamics 
simulations  are  carried  out  for  400  ps  and  the  results  used 
to  calculating  the  heat  current  for  every  time  step. 


3.  Results  and  discussion 

3.1.  The  effect  of  computational  cell  size 

As  discussed  earlier,  the  relative  magnitudes  of  the  com¬ 
putational  cell  size  with  respect  to  the  phonon  mean  free 
path  can  be  an  important  factor  affecting  the  accuracy  of 
thermal  conductivity  computed  using  atomistic  simulations. 
When  the  simulation  cell  is  too  small,  the  time  for  phonons 
to  travel  through  the  simulation  cell  is  shorter  than  the  de¬ 
cay  time  of  the  heat-current  auto-correlation  function.  This 
causes  phonon  scattering  to  take  place  more  frequently  than 
they  would  in  an  infinite  system.  In  such  cases,  only  the  short 
time  correlation  functions  are  expected  to  be  accurate.  Nev¬ 
ertheless,  Che  et  al.  [2]  showed  that  thermal  conductivity 
can  be  computed  using  computational  cells  smaller  than  the 
phonon  mean  free  path.  Using  the  macroscopic  laws  of  re¬ 
laxation  and  the  Onsager’ s  postulate  for  microscopic  thermal 
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Table  2 


Optical  and  acoustic  weighting  factors,  A0  and  Aa,  and  relaxation  times,  rG  and  ra,  obtained  by  nonlinear  least  squares  fitting  of  the  heat  current 
autocorrelation  function  for  (10,  10)  armchair  single-walled  carbon  nanotubes 


Numbers  of  atoms 

Cell  length  (nm) 

A0  (xlO~6J2/cm2) 

o 

43 

Aa  (xl0~6J2/cm2) 

ra  (ps) 

A0r0/A0  t„  +  Aara 

400 

2.477 

3.339 

0.074 

1.211 

17.597 

0.0012 

800 

4.954 

3.334 

0.070 

1.216 

17.870 

0.001 1 

1600 

9.908 

3.299 

0.072 

1.227 

17.940 

0.001 1 

3200 

19.816 

3.376 

0.069 

1.232 

17.851 

0.0011 

6400 

39.632 

3.339 

0.070 

1.230 

17.877 

0.0011 

Average 

3.343 

0.070 

1.229 

17.869 

fluctuations,  Che  et  al.  [2]  first  showed  that  the  long-time 
asymptotic  decay  of  the  heat-current  auto-correlation  func¬ 
tion  (which  primarily  controls  thermal  conductivity)  is  of  an 
exponential  type.  Consequently,  relatively  short  time  simu¬ 
lation  data  obtained  using  medium-size  computational  cells 
can  be  used  to  determine  quite  accurately  the  exponential  de¬ 
cay  parameters  of  the  heat-current  autocorrelation  function. 

To  determine  the  effect  of  the  computational  cell  size  on 
lattice  thermal  conductivity,  five  different  sizes  of  the  com¬ 
putational  cell  are  used  for  each  of  the  three  carbon  nan¬ 
otubes  analyzed  in  the  present  work.  The  number  of  atoms 
in  the  computational  cell  in  each  case  is  given  in  Tables  2-4. 
Due  to  the  fact  that  the  carbon  nanotubes  are  single  walled 
(one  atomic  layer  thick),  only  thermal  conductivity  in  the 
direction  of  nanotubes  axis  is  computed.  This  is  achieved  by 
using  only  the  data  for  thermal-energy  current  in  the  direc¬ 
tion  of  nanotubes  axis  when  calculating  the  auto-correlation 
functions,  Eq.  (5). 

Examples  of  typical  heat-current  auto-correlation  func¬ 
tions  for  the  (10,  10),  (18,  0)  and  (14,  6)  nanotubes  are 
shown  in  Fig.  2(a)-(c),  respectively.  These  function  are  ob¬ 
tained  using  computational  cells  containing  3200,  2880  and 


2528  atoms,  respectively.  In  general,  the  auto-correlation 
functions  are  characterized  by  a  rapid  initial  decay  followed 
by  a  gradual,  long-time  exponential  decay.  The  initial  fast 
decay  can  be  attributed  to  high-frequency  optical  phonon 
modes  which  are  associated  with  out-of-phase  vibrations  of 
the  atoms  residing  on  two  sub-lattices  in  the  nanotube  crystal 
structure.  These  phonons  are  scarcely  populated  and  weakly 
coupled  with  low-frequency  acoustic  (in-phase  vibrational) 
modes  at  room  temperature  and,  hence,  they  do  not  sig¬ 
nificantly  contribute  to  thermal  conductivity.  The  long-time 
behavior  of  the  heat-current  auto-correlation  functions  and, 
hence,  thermal  conductivity  are  controlled  by  low-frequency 
acoustic  phonon  modes. 

The  auto-correlation  function  results  such  as  the  one 
shown  in  Fig.  2(a)-(c)  are  fitted  using  the  Levenberg- 
Marquardt  nonlinear  least-squares  method  [23]  to  the  fol¬ 
lowing  double  exponential  function: 

Cj(t)  —  A0  exp  f  ^  exp  f  ^  ,  t  >  0,  (16) 

where  the  subscript  o  and  a  are  used  to  denote  the  opti¬ 
cal  and  acoustic  phonon  modes,  respectively.  Substitution 


Table  3 


Optical  and  acoustic  weighting  factors,  A0  and  Aa,  and  relaxation  times, 
autocorrelation  function  for  (18,  0)  zig-zag  single- walled  carbon  nanotubes 

To  and  ra, 

obtained  by  nonlinear  least 

squares 

fitting  of  the  heat  current 

Numbers  of  atoms 

Cell  length  (nm) 

A0  (xl0~6J2/cm2) 

r0  (ps) 

Aa  (xlO“6J2/cm2) 

ra  (ps) 

A0t0/A0t0  +  Aara 

360 

2.145 

3.327 

0.075 

1.135 

17.230 

0.0013 

720 

4.290 

3.329 

0.071 

1.141 

17.306 

0.0012 

1440 

8.580 

3.338 

0.073 

1.150 

17.395 

0.0012 

2880 

17.160 

3.330 

0.071 

1.153 

17.349 

0.0012 

5760 

34.320 

3.344 

0.069 

1.150 

17.401 

0.0012 

Average 

3.338 

0.070 

1.150 

17.375 

Table  4 

Optical  and  acoustic  weighting  factors,  A0 

and  Aa,  and  relaxation  times, 

x0  and  Ta, 

obtained  by  nonlinear  least 

squares 

fitting  of  the  heat  current 

autocorrelation  function  for  (14,  6)  single-walled  carbon  nanotubes 

Numbers  of  atoms 

Cell  length  (nm) 

A0  (xl0~6J2/cm2) 

To  (ps) 

Aa  (x  10-6  J2/cm2) 

ra  (ps) 

A0r0/A0r0  +  Aara 

632 

3.813 

3.280 

0.068 

1.129 

16.774 

0.0012 

1264 

7.626 

3.294 

0.070 

1.133 

16.803 

0.0012 

2528 

15.252 

3.289 

0.069 

1.134 

16.906 

0.0012 

5056 

30.504 

3.330 

0.065 

1.135 

16.899 

0.001 1 

Average 

3.311 

0.069 

1.131 

16.880 
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Fig.  2.  Heat-current  auto-correlation  functions  for:  (a)  a  (10,  10);  (b)  a  (18,  0)  and  (c)  a  (14,  6)  single-walled  carbon  nanotube  at  300  K. 


of  Eq.  (16)  in  Eq.  (4)  yields  the  following  expression  for 
thermal  conductivity: 

1 

1  =  kBT  2V('A°X0  +  Aata'>'  (17) 

Following  the  procedure  suggested  by  Che  et  al.  [2], 
the  auto-correlation  function  results  for  all  simulation  runs 
corresponding  to  the  first  3ps  of  the  correlation  time  are 
fitted  to  the  function  defined  in  Eq.  (16)  to  determine  the 
parameters  A0,  r0,  Aa,  and  ra.  The  results  of  this  proce¬ 
dure  are  given  in  Tables  2-4.  A  simple  analysis  of  the 
results  shown  in  Tables  2-4  indicates  that  the  contribution 
of  high-frequency  optical  phonon  modes  to  thermal  con¬ 


ductivity,  (A0t0)/(A0t0  +  Aara),  is  indeed  very  small  is 
typically  around  0.1%. 

The  dependence  of  thermal  conductivity  in  the  three  types 
of  carbon  nanotubes  analyzed  in  the  present  work  on  the  size 
of  the  simulation  cell  is  displayed  in  Fig.  3(a)-(c).  The  error 
bars  shown  in  Fig.  3(a)  and  (b)  correspond  to  ±1  standard 
deviation  for  the  results  of  five  molecular  dynamics  runs. 
The  results  displayed  in  Fig.  3(a)-(c)  show  that,  as  expected, 
when  the  simulation  cell  is  too  small  (i.e.  contains  less  than 
~1000  atoms),  the  atoms  in  a  region  of  the  simulation  cell 
do  not  have  enough  time  to  lose  their  previous  dynamic 
information  before  a  periodically  equivalent  phonon  arrives 
in  this  region.  Consequently,  since  the  corresponding  corre- 
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Fig.  3.  Lattice  thermal  conductivity  in:  (a)  a  (10,  10);  (b)  a  (18,  0)  and  (c)  a  (14,  6)  single-walled  carbon  nanotube  at  300  K. 


lation  functions  are  contaminated  by  such  memory  effects 
and  they  do  not  reflect  behavior  of  the  real  system,  com¬ 
puted  thermal  conductivity  is  not  very  accurate.  On  the  other 
hand,  when  the  computational  cell  size  is  sufficiently  large, 
computed  thermal  conductivity  is  accurate  and,  essentially, 
independent  of  the  size  of  the  computational  cell.  The  results 
shown  in  Fig.  3(a)-(c)  suggest  that  the  minimum  “critical” 
size  of  the  computational  cell,  corresponding  to  the  max¬ 
imum  correlation  time  of  3ps,  is  around  15-20  nm  and  is 
associated  with  cells  containing  approximately  3000-3500 
atoms. 


The  analysis  presented  above  established  a  minimum  crit¬ 
ical  computation  cell  size  beyond  which  the  cell  size  does 
not  affect  thermal  conductivity.  It  is  interesting  to  determine 
how  the  critical  cell  size  compares  with  the  phonon  mean 
free  path  in  single-walled  carbon  nanotubes.  To  estimate  the 
phonon  mean  free  path,  L,  the  following  relation  form  the 
kinetic  theory  of  solids  is  used  [20]: 

X  «  \CvpvL  (18) 

where  Cv  is  the  mass-based  constant-volume  specific  heat 
and  v  is  the  speed  of  sound.  The  molecular  dynamics 


M.  Grujicic  et  al.  / Materials  Science  and  Engineering  B107  (2004)  204-216 


211 


simulations  carried  out  in  the  present  work  yielded  the 
mean  interatomic  spacing  of  0.142nm.  The  density  of 
single-walled  carbon  nanotubes  is  difficult  to  determine 
since  they  are  one-atom  thick.  Since  single-walled  carbon 
nanotubes  are  typically  bundled  in  ropes  with  a  triangular 
arrangement  of  the  nanotubes  and  an  inter-tube  spacing 
equal  to  the  Van  der  Walls  radius  of  carbon  (0.17  nm),  the 
nanotubes’  wall  thickness  is  set  equal  to  this  value.  This  pro¬ 
cedure  yielded  the  density  of  p  =  2.3  g/cm3.  The  experimen¬ 
tal  values  for  the  specific  heat  of  Cv  =  500  J/(kg  K)  and  for 
compressibility  of  fj  —  0.024  GPa-1  for  single-walled  car¬ 
bon  nanotubes  of  comparable  diameters  has  been  taken  from 
Refs.  [21,22].  Using  the  following  equation:  v  =  1 
the  speed  of  sound  has  been  computed  as  4256  m/s,  Lastly, 
using  the  average  thermal  conductivity  for  the  three  types 
of  nanotubes,  X  =  16.5W/(cmK),  and  Eq.  (18),  the  mean 
free  path  for  acoustic  phonons  of  L  =  lOOOnm  has  been 
obtained.  Thus,  the  acoustic  phonon  mean  free  path  is  larger 
the  critical  computational  cell  size  by  a  factor  of  50-65. 

The  finding  presented  above  shows  that  despite  the  fact 
that  the  phonon  mean  free  path  is  considerably  larger  than 
the  sizes  of  computational  cells  used,  convergence  in  ther¬ 
mal  conductivity  can  be  obtained.  Fig.  3(a)-(c).  It  should  be 
noted  that  in  order  to  set  the  computational  cell  size  compa¬ 
rable  with  the  phonon  mean  free  path,  computational  cells 
containing  on  the  order  of  107— 108  atoms  would  have  to  be 
used.  Molecular  dynamics  simulations  using  such  a  large 
number  of  atoms,  while  feasible  particularly  if  an  advantage 
is  taken  of  the  parallel  computing,  would  be  computational 
very  expensive  and  are  not  very  appealing.  Instead,  as  ini¬ 
tially  shown  by  Che  et  al.  [2]  and  also  confirmed  in  the 
present  work,  smaller  size  simulation  cells  and  short-time 
heat-current  auto-correlation  functions  can  be  used  to  deter¬ 
mine  thermal  conductivity. 

The  values  of  the  weighting  factors  A0  and  ,4a  and  the  cor¬ 
responding  exponential  decay  constants  r0  and  ra  obtained 
by  the  nonlinear  least-squares  fitting  of  the  heat-current 
auto-correlation  functions  in  the  three  types  of  nanotubes 
and  for  the  simulation  cells  with  different  numbers  of  atoms 
are  listed  in  Tables  2-A.  To  determine  the  average  values  of 
these  parameters,  the  weighting  factors  A0  and  ,4a  are  nor¬ 
malized  by  the  nanotube  volume,  weighted  by  the  number  of 
atoms  in  each  simulation  and  averaged  over  all  simulations 
for  a  given  type  of  nanotube.  Similar  averaging  without  vol¬ 
ume  normalization  is  used  to  obtain  average  values  for  the 
relaxation  times,  r0  and  ra.  The  results  of  this  procedure  are 
listed  in  Tables  2—4  in  the  row  denoted  as  “average”. 

3.2.  The  effect  of  chirality 

The  results  displayed  in  Fig.  3(a)-(c)  further  show  that 
chirality  has  an  effect  on  the  lattice  part  of  thermal  con¬ 
ductivity  in  single-walled  carbon  nanotubes.  Specifically, 
thermal  conductivity  is  the  highest  {X  =  17.8W/(cmK)) 
in  the  (10,  10)  arm-chair  nanotube  and  the  lowest  (X  = 
15.6W/(cmK))  in  the  (14,  6)  nanotube  of  general  chiral¬ 


ity.  Since  these  results  pertain  to  the  lattice  contribution  to 
thermal  conductivity  alone,  they  have  to  be  attributed  to 
differences  in  phonon-phonon  interactions  and  the  resulting 
differences  in  mean  free  path  in  the  three  types  of  nanotubes. 
As  stated  earlier,  detailed  modeling  phonon-phonon  inter¬ 
actions  is  very  complicated  [13]  and  it  is  beyond  the  scope 
of  the  present  study.  Nevertheless,  it  is  well  established  that 
lattice  thermal  conductivity  is  controlled  by  the  so-called 
Umklapp  phonon  collisions  represented  as:  k\  +  kj  = 
k-i  +  G,  where  k\  and  ki  are  wave  vectors  of  the  colliding 
phonons,  kj,  the  wave  vector  of  the  resulting  phonon  and  G 
is  2tt  times  a  reciprocal  lattice  vector.  Due  to  differences 
in  chirality  and  the  resulting  differences  in  magnitudes  of 
the  periodic  length  in  the  axial  direction,  one  could  expect 
differences  in  the  permissible  G  vectors  and,  hence,  lattice 
thermal  conductivity  in  the  three  nanotubes  analyzed. 

It  should  be  noted  that,  based  on  the  magnitude  of  the 
band  gap  alone,  the  electronic  contribution  to  thermal  con¬ 
ductivity  can  also  be  expected  to  be  the  highest  in  the  (10, 
10)  nanotube  and  the  lowest  in  the  ( 14,  6)  nanotube.  It  should 
be  noted,  however,  that  in  addition  to  the  band  gap,  the 
electronic  band  stmcture  as  well  as  electron-electron  and 
electron-phonon  scattering  also  affect  the  electronic  contri¬ 
bution  to  thermal  conductivity.  Nevertheless,  the  observation 
that  lattice  thermal  conductivity  can  vary  by  as  much  as  20% 
with  nanotube  chirality  and  that  electronic  thermal  conduc¬ 
tivity  can  be  affected  in  the  same  direction  by  chirality,  sug¬ 
gests  that  thermal  conductivity  of  individual  single-walled 
nanotubes  in  nanotube  ropes  (consist  of  nanotubes  of  vari¬ 
ous  chirality)  can  vary  substantially  from  one  nanotube  to 
the  other. 

3.3.  The  effect  of  temperature 

The  molecular  dynamics  based  procedure  for  compu¬ 
tation  of  the  thermal  conductivity  in  single-walled  carbon 
nanotubes  described  in  Sections  2.1  and  2.2  is  utilized  in 
this  section  to  determine  the  temperature  dependence  of 
thermal  conductivity  in  these  materials.  Over  the  last  few 
years,  there  has  been  a  number  of  experimental  and  theoret¬ 
ical  investigations  of  the  effect  of  temperature  on  the  ther¬ 
mal  conductivity  of  single-walled  carbon  nanotube  bundles 
and  of  individual  multi-walled  carbon  nanotubes  [31-35]. 
However,  no  reliable  experimental  data  presently  exist  for 
the  thermal  conductivity  of  individual  single-walled  carbon 
nanotubes. 

The  effect  of  temperature  in  a  range  between  50  and 
400  K  on  the  mean  value  of  thermal  conductivity  in  the  three 
types  on  single-walled  carbon  nanotubes  analyzed  in  the 
present  work  is  shown  in  Fig.  4.  At  the  lowest  temperature 
explored,  thermal  conductivity  increases  with  an  increase  in 
temperature  while  in  the  upper  portion  of  the  temperature 
range  examined,  thermal  conductivity  decreases  with  an 
increase  in  temperature.  This  finding  suggests  that,  at  low 
temperatures,  thermal  conductivity  is  dominated  by  acous¬ 
tic  phonons  while,  at  high  temperatures,  phonon-phonon 
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Fig.  4.  The  effect  of  temperature  on  the  thermal  conductivity  in  the  three 
types  of  single-walled  carbon  nanotubes  analyzed  in  the  present  work. 


umklapp  scattering  causes  the  thermal  conductivity  to  de¬ 
crease  with  an  increase  in  temperature  [32]. 

For  comparison,  experimental  results  pertaining  to  the  ef¬ 
fect  of  temperature  on  the  thermal  conductivity  in  a  bundle 
of  the  single-walled  carbon  nanotubes  are  also  displayed  in 
Fig.  4  [34].  While  the  agreement  between  the  computational 
and  the  experimental  results  in  the  common  temperature 
range  can  be  characterized  as  only  fair,  the  true  validation 
of  the  present  computational  model  entails  the  experimen¬ 
tal  thermal  conductivities  of  isolated  single-walled  carbon 
nanotubes  with  the  known  chirality.  Unfortunately,  as  stated 
earlier,  such  experimental  data  are  presently  not  available. 

3.4.  A  comparison  with  the  Boltzmann  transport  equation 

In  this  section,  the  equation  of  phonon  radiative  transport 
(EPRT)  [29],  derived  from  the  Boltzmann  transport  equa¬ 
tion  [30]  is  used  to  compute  lattice  thermal  conductivity  of 
carbon  nanotubes  and  compare  it  with  the  results  presented 
in  the  previous  sections.  Since  details  of  this  calculation 
will  be  presented  in  a  future  communication,  only  a  brief 
overview  of  the  EPRT  approach  will  be  given  here.  For  a 
one -dimensional  case,  the  EPRT  equation  is  defined  as: 

13/®  ,  31a,  0.5/V®d/*-/® 

+  F— =  -  (19) 

v  at  ax  utr 

where  ho  is  the  flux  of  energy  in  the  direction  of  phonon 
propagation  per  unit  area,  per  unit  solid  angle,  per  unit  fre¬ 
quency,  t  is  time,  x  the  nanotube  axial  direction,  co  frequency 
and  p.  —  cos  9  (9  the  polar  angle),  tr  the  phonon  relaxation 
time. 

Eq.  (19)  is  applied  to  the  case  of  a  10  pm  long  nanotube, 
subjected  at  t  —  0  to  a  temperature  difference  of  0. 1  K  be¬ 
tween  its  ends  and  solved  using  a  first-order  upward  finite 


Normalized  Distance  Along  Nanotube  Length 


Normalized  Distance  Along  Nanotubc  Length 

Fig.  5.  Steady-state  temperature  (a)  and  heat  flux  (b)  distributions  along 
the  nanotube  length  computed  using  the  Equation  of  Phonon  Radiative 
Transfer  (29). 


difference  method.  In  these  calculation,  tr  =  ra  =  17  ps  is 
used.  The  calculations  of  I  a,  are  carried  out  until  a  steady 
state  is  reached  and  then  to  integrations  procedures  are  ap¬ 
plied  to  compute  the  temperature  and  the  heat  flux.  The  re¬ 
sulting  steady  state  temperature  and  heat  flux  distributions 
along  the  nanotube  length  are  shown  in  Fig.  5(a)  and  (b),  re¬ 
spectively.  It  should  be  noted  that  a  normalized  temperature 
axis,  (7(K)  —  300.0) /0.1,  is  used  in  Fig.  5(a)  and  normalized 
length  axes,  x  (pm)/0.1,  are  used  in  Fig.  5(a)  and  (b). 

Using  the  results  shown  in  Fig.  5(a)  and  (b)  and  a 
one -dimensional  form  of  the  steady-state  Fourier  heat- 
conduction  equation,  Eq.  (1),  thermal  conductivity  is 
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computed  as  8.3  W/(cmK).  This  value  is  smaller  by  a  factor  Within  the  AIREBO  potential  formulation,  the  potential 

of  ~2  than  its  counterparts  displayed  in  Fig.  3(a)-(c).  energy  of  a  system  of  atoms  is  represented  as  [17]: 


4.  Conclusions 

Based  on  the  results  obtained  in  the  present  work,  the 

following  main  conclusions  can  be  drawn: 

1.  Due  to  an  exponential-decay  character  of  the  long-time 
current-energy  auto-correlation  functions,  the  lattice 
contribution  to  thermal  conductivity  of  single-walled 
carbon  nanotubes  can  be  determined  quite  accurately  us¬ 
ing  molecular  dynamics  simulations  and  computational 
cells  substantially  smaller  than  the  phonon  mean  free 
path. 

2.  The  lattice  contribution  to  thermal  conductivity  can  vary 
by  as  much  as  20%  in  single-walled  carbon  nanotubes 
depending  on  their  chirality. 

3.  The  lattice  contributions  to  thermal  conductivity  com¬ 
puted  using  molecular  dynamics  based  thermal-current 
auto-correlation  functions  and  the  ones  computed  using 
the  Boltzmann  transport  equation  agree  within  a  factor 
of  two  with  each  other. 
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Appendix  A.  The  adaptive  intermolecular  reactive 
empirical  bond  order  (AIREBO)  potential 

The  Brenner’s  Reactive  Empirical  Bond  Order  (REBO) 
potential  [18,24-26]  is  an  interactive  potential  initially  de¬ 
veloped  to  model  the  covalent  bonding  interactions  in  carbon 
and  hydrocarbon  systems  attending  Chemical  Vapor  Depo¬ 
sitions  (CVD)  of  diamond.  While  this  potential  and  its  vari¬ 
ous  incarnations  have  been  recently  extended  to  analyze  en¬ 
ergetic,  elastic  and  vibrational  properties  of  fullerenes  [26], 
carbon  nanotubes  [27],  amorphous  carbon  [28],  etc.  they 
are  not  the  most  appropriate  for  such  analyses  due  to  sig¬ 
nificant  contribution  of  nonbonding  (intermolecular)  inter¬ 
actions  which  the  REBO  potential  does  not  include.  Re¬ 
cently,  Stuart  et  al.  [17]  upgraded  the  REBO  potential  by 
including  an  adaptive  treatment  of  the  interactions  involv¬ 
ing  nonbonded  atoms  and  dihedral-angle  intermolecular  in¬ 
teractions.  This  potential,  named  the  Adaptive  Intermolecu¬ 
lar  Reactive  Empirical  Bond  Order  (AIREBO)  potential  is 
briefly  overviewed  in  this  section. 


‘  j¥=i 


Z7  REBO 
Eij 


+  dGE  E  £“ 

k^i,j 


(A.l) 


The  REBO  potential,  ,  describes  the  interactions  be- 

tween  covalently-bonded  atoms  (the  intramolecular  inter¬ 
actions),  while  the  dispersion  and  intermolecular  repulsion 
interactions  between  nonbonded  atoms  are  accounted  for 
through  the  use  of  the  Lennard-Jones  (LJ)  term.  The  last 
term  in  Eq.  (A.l)  is  used  to  describe  dihedral-angle  inter¬ 
molecular  interactions  which  are  deemed  not  significant  in 
the  analysis  of  carbon  nanotubes  [17]  and,  hence,  this  term 
is  not  considered  in  the  present  work. 

The  REBO  interaction  term,  £EEB0,  is  defined  as  [18]: 

Ef] 80  =  If  +  bijVfr  (A.2) 


where  the  repulsive,  V?,  and  the  attractive,  V/\  interaction 
energies  between  atoms  i  and  j  separated  by  r,;  in  Eq.  (A.2) 
are  combined  in  a  ratio  determined  by  the  many-body  bond¬ 
ing  parameter,  bij.  The  repulsive  term  is  defined  as  [18]: 


vij  =  wu(nj) 


l  +  Ql 

rij  . 


Ajj  e  “Brw 


(A. 3) 


where  the  parameters  Qjj,  Ajj,  and  a,j  depend  on  the  atom 
types  i  and  j.  Values  for  these  and  all  other  REBO  potential 
parameters  can  be  found  in  Table  2,  Ref.  [17].  The  w\]  term 
in  Eq.  (A. 3)  is  a  bond-weighting  function  defined  as: 

Wjjinj)  =  S'Mnj)),  (A.4) 


which  is  used  to  switch  off  the  REBO  interactions,  in  a 
continuous  manner,  when  the  distance  between  atoms  i  and 
j  exceeds  an  atomic-type  dependent  bonding  distance  and 
the  i-j  bond  breaks.  The  switching  function,  S'(tc(rjj)),  and 
the  scaling  function,  tc(rjj),  are  defined  as: 

S' {t)  =  ©(-f)  +  ©(f) 0(1  -  t)\[  1  +  cos  (;rr)]  (A.5) 


and 

.  _  7-™° 

.  ,  riJ  rij 

tc(rjj)  =  - ; — - 

J  max  min 

ij  ij 


(A. 6) 


respectively,  while  the  switching  region  is  defined  in  terms 
of  r™m  and  r™ax.  ©  in  Eq.  (A.5)  denotes  the  Heaviside  step 
functions.  S' (t)  is  unity  for  t  <  0  and  zero  for  t  >  1  and 
provides  a  smooth  switching  between  these  two  values  for 
0  <  t  <  1. 

The  attractive  pair  interaction  term  in  Eq.  (A.2)  is  given 
by  a  triple  exponential  function  as  [18]: 

V*  =  (A.7) 

n=  I 


and  is  switched  off  smoothly  over  the  same  switching  re¬ 
gion  as  Vj)  through  the  use  of  the  bond  weighting  function 
WiAnj). 
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The  many-body  bonding  term,  bjj,  in  Eq.  (A. 2)  defines  the 
“bond  order”  for  the  interaction  between  atoms  i  and  j  (the 
larger  bjj,  the  stronger  the  i-j  bond),  by  making  the  strength 
of  the  i-j  bond  dependent  on  the  local  atomic  environment 
and  is  given  as: 

bij  =  \\pp  +  PjH  +  rfj  +  •'A-8) 

The  principal  contribution  to  bij  in  Eq.  (A. 8)  comes  from 
the  covalent-bond  interaction  terms,  p™ ,  and  pjj1  (p™  is 
not  necessarily  equal  to  Pp),  where  p™  is  defined  as: 

-  -1/2 

1  +  £  Wij(rik)gi(cosOjik) eXjik  +  Py  (A.9) 
k&J 

where  Qjjk  is  the  bond  angle  between  the  rj,  vector  connecting 
neighboring  atoms  j  and  i  and  the  vectors  /y,,-  connecting 
atom  i  to  another  neighboring  atom  k.  The  penalty  function, 
gi,  in  Eq.  (A.9),  is  used  to  provide  an  increase  in  the  energy 
for  bonds  that  are  too  close  to  one  another  and  its  functional 
form  is  a  fifth-order  spline.  The  coefficients  of  the  spline  are 
listed  in  Table  7,  Ref.  [17]. 

The  two  remaining  terms  in  Eq.  (A.9)  represent  relatively 
small  correction  factors.  The  ex’ik  term  is  added  to  improve 
the  potential  energy  surface  for  abstraction  of  hydrogen 
atoms  from  hydrocarbons  while  the  Pjj  term  (a  two  dimen¬ 
sional  cubic  spline)  is  used  to  obtain  accurate  bond  energies 
for  small  hydrocarbons.  The  parameters  for  these  two  terms 
are  given,  respectively,  in  Tables  2  and  8,  Ref.  [17]. 

In  addition  to  the  covalent-bonding  interactions  given  by 
Eq.  (A.9),  the  REBO  potential  also  includes  contributions  to 
the  bond  order  from  radical  and  conjugation  effects.  These 
enter  the  potential  through  the  jv'k  term  in  Eq.  (A.  8),  which 

is  a  three  dimensional  cubic  spline  of  Nij,  Njj,  and  (V™nj  as 
independent  variables.  Nij  and  /V;,  are,  respectively,  atoms  i 
and  j  coordination  numbers,  while  A™nj  is  a  local  measure 
of  conjugation  in  the  i-j  bond  and  is  given  as: 


^Conj  _  j 


-i  2 


£  SkcWik(nk)S'  (^conj  (Nki)) 
k^ij 

2 

£  8lcWjl(rjl)S'  (/conj  (Nij)) 
¥>J 


(A.  10) 


fconj  in  Eq.  (A.  10)  specifies  the  range  of  coordination  num¬ 
bers  under  which  a  bond  is  assumed  to  be  part  of  a  radical 
or  conjugated  network  and  is  defined  as: 


entirely  on  local  coordination.  The  interpolation  points  for 
the  three-dimensional  spline,  tTc,  are  provided  in  Table  9, 
Ref.  [17], 

The  last  contribution  to  the  bond  order  bij  in  Eq.  (A. 8)  is 
7 rjjh .  This  term  imposes  a  penalty  for  rotation  around  multiple 
bonds  and  is  defined  as: 


-  Tii(Nih  Nju  Np)  £  £(1-  cos2  cokiji) 

k^iJljiiJ 

x  w'ik(rikWji(rjl)0(sin(9jik)  -  smm ) 
x  0(sin(0y/)  -  .vmin)  (A.  12) 


where  Tjj  is  a  three-dimensional  cubic  spline,  with  interpo¬ 
lation  points  given  in  Table  10,  Ref.  [17].  The  torsion  angle 
cokiji  is  defined  as  the  angle  between  the  plane  defined  by  the 
vectors  /-,/■  and  r(/  and  that  defined  by  the  vectors  r,;  and  rji 
as: 


COS  U>kijl  — 


rji  x  rik 
I  rji  x  rik  | 


nj  x  rH 

\rij  x  rji\ 


(A.  13) 


The  bond  weighting  function,  w'iry),  used  in  Eq.  (A.  12) 
is  defined  as: 


wUnj)  =  S'(t'c{rij)) 


(A.  14) 


and  differs  slightly  from  that  defined  in  Eq.  (A. 4),  through 
a  different  scaling  function  t'c  defined  as: 


t'c  (nj)  = 


min 

ij 


rmax 

ij 


(A.  15) 


The  El.j]  term  in  Eq.  (A.l)  is  based  on  the  Lennard- Jones 


(LJ)  12-6  potential  and  is  defined  as: 


Vipnj)  =  4eij 


12 


(A.  16) 


The  inclusion  of  the  LJ  term  creates  several  challenges, 
foremost  of  these  is  introduction  of  a  steep  repulsive  wall, 
which  prevents  unbonded  atoms  (atoms  associated  with  dif¬ 
ferent  molecules)  from  getting  close  enough  to  interact  via 
the  REBO  potential.  Rather  than  overcoming  this  problem 
by  incorprating  a  simple  distance-dependent  switching  func¬ 
tion  in  the  LJ  potential,  and  thus  neglecting  the  effect  of 
chemical  environment,  Stuart  et  al.  [17]  introduced  three  cri¬ 
teria  to  switch  on  or  off  the  LJ  interactions  in  an  adaptive 
fashion.  The  three  criteria  are  based  on:  (a)  the  separation 
distance  between  the  two  atoms  in  question;  (b)  their  bond 
order;  and  (c)  the  network  of  covalent  bonds  connecting  the 
atoms.  Consequently  the  LJ  interaction  term  in  Eq.  (A.l)  is 
defined  as: 


Wcoiy  _  „conj,min 

fconj  (  Np )  =  - j - ^ (A.ll) 

J  lJ  ,rconi,max  ,rconi,min  v  7 

Nij  -  Nij 

The  Np  variable,  which  is  unity  for  nonconjugated  bonds 
and  can  be  as  high  as  nine  in  polyaromatic  compounds, 
is  an  empirical  measure  of  unsaturation  that  is  dependent 


4J  =  SUAnjjiSiMb^CyV^inj) 

+  [1  -  S{tt(fij))]CijV]j\rij),  (A.  17) 

where  the  atomic-separation-distance  based  switching  func¬ 
tion,  S(t),  is  given  as: 

S(t )  =  6>(-f)  +  <9(f)0(  1  -  r)[l  -  f2(3  -  2 r)]  (A.  18) 
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and  the  scaling  function,  tT(rij),  in  Eq.  (A. 17)  is  given  as: 


„  ,,LJ  min 

r,J  rij 

LJ  max  LJ  min 

ij  ij 


(A.  19) 


When  rij  >  fTJmax,  5(fr(r,y))  =  0  and  the  magnitude  of 
the  i-j  atomic  separation  distance,  r(/,  has  no  effect  on  the 
LJ  interactions.  When  r,,  <  r|;J  max ,  on  the  other  hand,  the 
first-term  in  Eq.  (A.  17)  is  non  zero  and  the  LJ  interactions 
are  contingent  on  the  values  of  S(t^(b*))  and  Cjj. 

The  switching  region  |>kJmin,  rUmaxj  jn  gq  (A. 19)  is 
chosen  such  that  no  artificial  reaction  barrier  is  introduced 
due  to  LJ  repulsions  at  short  distances  (rjjJ  "lln  =  cry)  and 
that  switching  off  leaves  the  minimum  of  the  LJ  potential 
function  unperturbed  (,TJmax  =  21/6tTy). 

The  second  (bond  strength)  criterion  affecting  the  LJ  inter¬ 
actions  is  included  through  the  use  of  the  S(t\,(b*-))  switch¬ 
ing  function  in  Eq.  (A.  17)  where  the  scaling  function  t^(bjj) 
is  given  as: 

l)*  —  /jI-J  min 

th(bij)  =  (A. 20) 

J  ^max  —  b™m 


When  the  two  atoms  in  question  are  covalently  bonded,  (the 
bond-order  parameter,  /?*,  is  large  and,  hence,  t\,  >  1  and 
S(tb(b*))  =  0),  the  repulsive  LJ  interactions  (represented  by 
the  first  term  in  Eq.  (A.  17))  vanishes.  Lor  lower  values  of  the 
bond  order  parameter,  on  the  other  hand,  the  LJ  interactions 
will  be  included  to  a  variable  degree  depending  on  the  value 
of  b*:.  Lor  sufficiently  small  values  of  b*:  (which  indicates 
that  covalent  bonding  is  not  likely),  S(tb(b*j))  =  1,  and  the 
LJ  repulsion  interactions  are  undiminished. 

It  should  be  noted  that  separations  between  nonbonded 
atoms  generally  exceed  the  maximum  covalent  bonding  dis¬ 
tance,  rj’lax  in  Eq.  (A. 6),  and  hence  the  REBO  parameter, 
bjj,  does  not  provide  an  accurate  representations  of  the  bond 
order.  To  overcome  this  problem,  the  AIREBO  bond  order 
parameter  for  non-bonded  i-j  atomic  interactions,  b*j,  is  rep¬ 
resented  by  a  hypothetical  REBO  bjj  term  evaluated  at  r'™n 
with  the  distances  between  atoms  i  and  j  to  each  of  their 
neighbors  being  kept  unchanged. 

The  third  criterion  governing  the  LJ  interactions  is  repre¬ 
sented  by  the  connectivity  switch  Cjj  in  Eq.  (A.  17)  and  is 
introduced  to  suppress  LJ  interactions  between  first  neigh¬ 
bors  (1-2),  second  neighbors  (1-3)  and  third  neighbors  ( 1 — 4) 
which  are  either  well  described  by  the  REBO  potential  or  by 
the  dihedral-angle  potential  described  below.  To  derive  an 
expression  for  Cjj,  a  bond  weight  function  Wy(r)y  is  first  de¬ 
fined  to  enable  a  smooth  transition  between  bonded  (wij  = 
1)  and  nonbonded  ( wy  =  0)  atomic  interactions.  The  bond 
weight  function  is  given  as: 

Wijirij)  =  S\tc(rij ))  (A. 21) 

where  the  switching  function  S\t )  is  given  as: 

S'(t)  =  0(-t)  +  0(f)0(l  -  i)\[\  +  cos  (itt)]  (A.22) 


while  the  scaling  function  fc(^y)  is  defined  as: 

..  ..min 
.  .  1  'J  rij 

tcifij)  =  - - — — 

v  J'  max  _  min 

ij  ij 


(A. 23) 


To  exclude  (1-2),  (1-3)  and  (1^1)  LJ  interactions,  the 
connectivity  switch  Cjj  is  defined  as: 


Cij  =  1  -  max{wy(ry),  wik{rik)wkj{rkj),  Vk; 

Wik(rik)wki(rki)wij(rij),  Vk,  1}  (A.24) 

Thus,  when  atoms  i  and  j  are  neighbors  (ivyiry)  —  1)  or 
connected  by  one  ( u’ik(rik)wkj(rkj )  =  1)  or  two  ( Wik(jik )  • 
Wkl(rkl)  ■  wij(rlj )  =  1)  neighbors,  Cy  =  0  and  there  are  no 
LJ  interactions  between  the  atoms. 

In  summary,  the  strength  of  LJ  interactions  is  affected 
by  atomic  separation  distance,  bond  order  and  connectivity. 
Lor  the  LJ  interactions  to  be  fully  included,  the  atom  pairs 
in  question  must  not  be  first,  second  or  third  neighbors  and 
must  either  be  beyond  a  cut-off  distance  or  have  a  very  low 
value  of  the  bond  order. 
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